%% 写这个是因为一次性可能跑不完，每次都重跑又太浪费了
clear
clc
close all
%% 基本设置
file = 'danji_fuzzy'; %simulink模块名
load_system([file,'.slx']);%加载simulink模块
fuzzymzh = readfis('fuzzymzh.fis');%读取模糊控制器
open_system([file,'.slx']);%打开simulink模型

%定义初始论域
Range.e = [-25 25];
Range.ec = [-15 15];
kp = 0.9014;alpha_p = 0.05;
ki = 0.0921;alpha_i = 0.2;
kd = 0.6991;alpha_d = 0.2;

Range.Kp = kp.*(1+alpha_p.*[-1 1]);%给个初始值
Range.Ki = ki.*(1+alpha_i.*[-1 1]);
Range.Kd = kd.*(1+alpha_d.*[-1 1]);

%定义论域比例的上下限
lb = [0.8,0.08,0.6,0.03,0.1,0.1];%下界
ub = [1.0,0.10,0.8,0.07,0.3,0.3];%上界

npop=100;%种群个数
nvar=6;%参数量7
maxit=200;%迭代次数
pc1=0.8;%交叉率
pc2=0.6;
pm1=0.2;%变异率
pm2=0.1;

rng(0);%固定随机数

template.x=[];
template.y=[];
template.dominationset=[];%支配集，指个体能够支配的其他个体所对应的下标
template.dominated=[];%被支配数，指个体在自然选择中被支配的次数
template.rank=[];%等级，指个体在自然选择中所占据的生态位
template.cd=[];%拥挤度
 
load pop100.mat;%读进来
%% 初始化模糊控制器
fuzzymzh.Inputs(1).Range = Range.e;
fuzzymzh.Inputs(2).Range = Range.ec;
e = linspace(Range.e(1), Range.e(2), 7);
ec = linspace(Range.ec(1), Range.ec(2), 7);
for i = 1:7
    fuzzymzh.Inputs(1).mf(i).Parameters = e(i) + [-diff(Range.e) / 6 0 diff(Range.e) / 6];
    fuzzymzh.Inputs(2).mf(i).Parameters = ec(i) + [-diff(Range.ec) / 6 0 diff(Range.ec) / 6];
end

%% 交叉变异迭代
figure('Name','遗传算法迭代过程')%迭代过程

y1_mean = [];%迭代过程记录
y2_mean = [];
iteration = [];

n_best = 0;
y_best = 0;

while n_best <= maxit/4 & length(iteration) <= maxit
    tic
    it = length(iteration)+1;
    npc=1;
    popc=repmat(template,npop/2,2);
    fall = [pop.y];%fall存储所有的y值
    fall = reshape(fall,[2, npop]);%第一行存储所有的y1，第二行存储所有的y2
    fmax =2; %这一代中综合y值最高的个体，以其作为归一化的上限
    fave =mean(fall(1,:))./max(fall(1,:))+mean(fall(2,:))./max(fall(2,:)); %这一代中综合y值的平均数
    for i=1:npop/2%交叉
        ind=randperm(npop,2);%选择交叉个体

        fcross1=pop(ind(1)).y(1)./max(fall(1,:))+pop(ind(1)).y(2)./max(fall(2,:));
        fcross2=pop(ind(2)).y(1)./max(fall(1,:))+pop(ind(2)).y(2)./max(fall(2,:));
        fcross = max([fcross1 fcross2]);
        pc=pc2;
        if(fcross>fave)
            pc = pc1-(pc1-pc2).*(fcross-fave)./(fmax-fave);
        end

        value = rand();
        if(value<=pc)
        [popc(npc,1).x,popc(npc,2).x]=Cross(pop(ind(1)).x,pop(ind(2)).x);%交叉得到新的参数组合
        popc(npc,1).y=respond(file,popc(npc,1).x,fuzzymzh,Range);  
        popc(npc,2).y=respond(file,popc(npc,2).x,fuzzymzh,Range);
        npc=npc+1;
        end
    end

    npc = npc-1;
    popc(npc+1:npop/2,:)=[];
    npm=1;
    popm=repmat(template,npop,1);
    i_mutate = [];%存储变异个体的编号
    for j=1:npop%
        ind=randperm(npop,1);%选择变异个体
        fmutate=pop(ind(1)).y(1)./max(fall(1,:))+pop(ind(1)).y(2)./max(fall(2,:));
        pm=pm2;
        if(fmutate>fave)
            pm = pm1-(pm1-pm2).*(fmutate-fave)./(fmax-fave);
        end
        value=rand();
        if(value<=pm)
            i_mutate = [i_mutate ind];%存储变异个体
        end
    end
    mutate_rand = lhsdesign(length(i_mutate),1);%根据变异个体数量进行拉丁超立方抽样
    for j = i_mutate
        popm(npm,1).x=Mutate(pop(j).x,lb,ub,mutate_rand(npm));
        popm(npm,1).y=respond(file,popm(npm).x,fuzzymzh,Range);
        npm=npm+1;
    end
    npm=npm-1;
    popm(npm+1:npop)=[];
    popc=popc(:);

    newpop=[pop;popc;popm];%新种群
    [newpop,F]=Non_dominate_sort(newpop);
    newpop=Crowd(newpop,F);%计算拥挤度
    newpop=nsga2Sort(newpop);%根据生态等级与拥挤度进行排序
    pop=newpop(1:npop);%只取前npop个个体存活至下一次迭代

    y1_mean = [y1_mean,mean(fall(1,:))];
    y2_mean = [y2_mean,mean(fall(2,:))];
    iteration = [iteration,it];

    subplot(3,1,1)
    plot(1./fall(1,:),1./fall(2,:),'r*');%y1是稳态误差，y2是调节时间
    numtitle=num2str(it);
    title(['迭代次数=',numtitle]);
    xlabel('y1');
    ylabel('y2');
    subplot(3,1,2)
    plot(iteration,1./y1_mean);
    subplot(3,1,3)
    plot(iteration,1./y2_mean);
    set(gcf,'color','white');
    pause(0.001);

    frame = getframe(gcf);
    I=frame2im(frame);
    [I,map]=rgb2ind(I,256);
    if it == 1
        imwrite(I,map,'test.gif','gif', 'Loopcount',inf,'DelayTime',0.08);
    else
        imwrite(I,map,'test.gif','gif','WriteMode','append','DelayTime',0.08);
    end
    y_sum = 1*(0.0249.*fall(1,:))+1*(17.6135.*fall(2,:));%括号外是权重，括号内是普通pid的性能
    y_best(it+1) = max(y_sum);
    if y_best(it+1) == y_best(it)
        n_best = n_best+1;
    else
        n_best = 0;%最优解不一样，则最优解重复次数归0
    end
    toc
    save pop100_1.mat pop;%把最后的种群保存下来
end
save pop100_1.mat pop;%把最后的种群保存下来
disp(['第',num2str(find(y_sum==y_best(it))),...
            '个个体最优，y_sum=',num2str(y_best(it))]);